A Pilot Study on Patient-specific Computational Forecasting of Prostate Cancer Growth during Active Surveillance Using an Imaging-informed Biomechanistic Model

Abstract Active surveillance (AS) is a suitable management option for newly diagnosed prostate cancer, which usually presents low to intermediate clinical risk. Patients enrolled in AS have their tumor monitored via longitudinal multiparametric MRI (mpMRI), PSA tests, and biopsies. Hence, treatment is prescribed when these tests identify progression to higher-risk prostate cancer. However, current AS protocols rely on detecting tumor progression through direct observation according to population-based monitoring strategies. This approach limits the design of patient-specific AS plans and may delay the detection of tumor progression. Here, we present a pilot study to address these issues by leveraging personalized computational predictions of prostate cancer growth. Our forecasts are obtained with a spatiotemporal biomechanistic model informed by patient-specific longitudinal mpMRI data (T2-weighted MRI and apparent diffusion coefficient maps from diffusion-weighted MRI). Our results show that our technology can represent and forecast the global tumor burden for individual patients, achieving concordance correlation coefficients from 0.93 to 0.99 across our cohort (n = 7). In addition, we identify a model-based biomarker of higher-risk prostate cancer: the mean proliferation activity of the tumor (P = 0.041). Using logistic regression, we construct a prostate cancer risk classifier based on this biomarker that achieves an area under the ROC curve of 0.83. We further show that coupling our tumor forecasts with this prostate cancer risk classifier enables the early identification of prostate cancer progression to higher-risk disease by more than 1 year. Thus, we posit that our predictive technology constitutes a promising clinical decision-making tool to design personalized AS plans for patients with prostate cancer. Significance: Personalization of a biomechanistic model of prostate cancer with mpMRI data enables the prediction of tumor progression, thereby showing promise to guide clinical decision-making during AS for each individual patient.


Introduction
Prostate cancer is the second most common type of cancer and the fifth leading cause of cancer death among men worldwide (1).The clinical management of prostate cancer relies on two key strategies (2): regular screening and patient triage into risk groups.Regular screening of prostate cancer is performed in men over 50 and consists of a digital rectal exam and the measurement of the serum level of the PSA.If the results of these two tests are suggestive of marker of prostate cancer aggressiveness and is defined upon the cancerous tissue architecture observed in biopsy samples or surgically-removed prostates (2,4).Thanks to the current clinical management of prostate cancer, most tumors are detected and successfully treated at early organ-confined stage, which usually poses a low to intermediate risk to the patient (5,6).Although standard treatments for prostate cancer (e.g., surgery, radiotherapy) exhibit high rates of overall and disease-free survival (2,7,8), many newly diagnosed prostate cancer cases are known to be indolent and may not produce any symptoms or require treatment for a long time.These patients are prone to potential treatment side effects (e.g., incontinence, impotence) that can adversely impact their quality of life without improving longevity.Thus, overtreatment of prostate cancer is a major concern in the current clinical management of the disease (3,9,10).In addition, undertreatment of prostate cancer due to limitations in diagnosis and management (e.g., biopsy sampling errors, mpMRI confounders, limited age and comorbidity assessment) constitutes another important clinical challenge, because it may result in rapid growth of aggressive tumors, treatment failure, and reduced survival (3,(9)(10)(11).
As an alternative to direct treatment after diagnosis, active surveillance (AS) is another standard-of-care management option that is suitable for many newly diagnosed prostate cancer cases exhibiting low to intermediate clinical risk (2,3,8,12).In AS, patients are closely monitored via longitudinal mpMRI scans, digital rectal exams, PSA tests, and biopsies, such that treatment is delayed until these monitoring tests confirm prostate cancer progression to higher risk.Thus, AS is widely regarded as an ideal clinical strategy to reduce the overtreatment of newly diagnosed prostate cancer (2,3,12).In addition, the changes in longitudinal data collected during AS (e.g., mpMRI-observable features, PSA kinetics) can be leveraged to inform clinical decision-making and, hence, contribute to reduce treatment excesses and deficiencies (12,13).However, current AS protocols largely rely on recommending monitoring tests either at fixed times (e.g., mpMRI scan potentially followed by biopsy every 6-36 months, PSA test every 3-12 months) or after the observation of clinical events suggesting tumor progression (e.g., imaging after a sustained PSA rise, biopsy after radiological worsening of an mpMRI-detected lesion; see refs.2,12,13).This population-based, observational approach largely ignores the complex underlying tumor dynamics and the intrinsic heterogeneity of prostate cancer between and within patients (3,4,(13)(14)(15), complicates the design of personalized plans to ensure optimally controlled tumor monitoring, and hinders the early detection of prostate cancer progression, which may occur between two consecutive and widely-spaced monitoring tests.Thus, the optimal timing and type of monitoring tests for each individual prostate cancer case constitutes a central challenge in AS, along with the definition of triggering criteria to switch from monitoring to radical treatment (2,3,12,13,16).To address these unmet and timely demands in AS, we propose to leverage personalized computational predictions of prostate cancer growth (17,18).
Computational tumor forecasting is a technology that enables the estimation of the growth of tumors and their response to treatment on a personalized basis (18).These tumor forecasts are obtained by leveraging mathematical models based on differential equations or agent-based formulations that describe the main biological mechanisms underlying the development of these diseases and the effect of specific treatments (18,19).By calibrating the model parameters governing tumor dynamics with patient-specific data, the ensuing personalized model enables the representation of observed tumor growth as well as the prediction of the future development of the tumor, which can be leveraged to guide monitoring strategies or adjust treatment regimens for each individ-ual patient (18)(19)(20)(21)(22)(23)(24)(25)(26).Indeed, the increasing success of computational tumor forecasting has motivated its use to design digital twin technologies to systematically guide clinical decision-making using a predictive, patient-specific approach (27).In particular, medical imaging data, such as mpMRI scans, provide spatially-resolved anatomic and physiologic information about tumors that enable the personalization of spatiotemporal biomechanistic models based on partial differential equations, which can then render three-dimensional (3D) tumor forecasts within the patient's affected organ (18).While this type of models has been extensively used to predict the development and therapeutic response of other tumors (20)(21)(22)(23)(28)(29)(30)(31), there is a dearth of applications for prostate cancer (17,(32)(33)(34)(35). Instead, most biomechanistic models developed for prostate cancer consist of time-resolved differential formulations primarily based on serum PSA dynamics that have been successfully used to predict response to specific treatments, especially in the advanced stage of the disease (24)(25)(26)36).The abundance of these PSA models is primarily motivated by the extensive use of this biomarker to monitor prostate cancer, while longitudinal mpMRI data are almost only available in AS.However, time-resolved models lack spatial information on prostate cancer growth that can be key during AS, for example, to plan biopsies or to assess malignancy in prostate areas that may associate with worse prognosis (e.g., prostate capsule, regions adjacent to neurovascular bundles and seminal vesicles; see refs.[2][3][4]13).
Here, we present a pilot study of personalized forecasting of untreated prostate cancer growth during AS using a spatiotemporal biomechanistic model informed by longitudinal mpMRI data collected for each individual patient.
Toward this end, we develop a clinical-computational pipeline to segment the prostate and tumor, perform intrascan and interscan registration of imaging data, build a 3D virtual representation of the prostate geometry holding the longitudinal spatial changes of the tumor, calibrate a biomechanistic model of prostate cancer growth, and perform patient-specific tumor forecasting (see Fig. 1).First, we analyze whether our computational framework can represent prostate cancer growth by informing the model with three imaging datasets for each patient.Then, we investigate whether we can predict prostate cancer growth after informing the model with two imaging datasets by comparing a prediction of the personalized model against a third mpMRI dataset.In addition, because prostate cancer risk and AS eligibility heavily rely on GS (2,12), we further analyze our personalized computer simulations to search for potential model-based predictive biomarkers of prostate cancer progression to higher-risk disease.

Patient Data
Anonymized patient data were retrospectively collected at The University of Texas Health Science Center San Antonio (San Antonio, TX) following an Institutional Review Board (IRB)-approved and Health Insurance Portability and Accountability Act-compliant protocol, which did not require informed consent.The inclusion criteria were: (i) newly diagnosed, untreated, organconfined prostate cancer managed with AS after diagnosis; (ii) availability of three mpMRI scans during AS; (iii) a visible tumor across all mpMRI scans; and (iv) histopathologically-confirmed prostate cancer via biopsy and/or analysis of the surgically-removed prostate after AS.The cohort leveraged in this study consisted of 7 patients (n = 7), each of them having one organ-confined tumor matching the inclusion criteria.The longitudinal mpMRI scans of each patient were collected over a period of median (range) of 3.4 (1.4,4.9) years with

FIGURE 1
Imaging-informed computational pipeline for tumor forecasting to support clinical decision-making in AS for prostate cancer.
A, Standard-of-care AS protocol for an illustrative patient.Following the detection of a serum PSA level moderately larger than 4 ng/mL, the patient undergoes an initial mpMRI scan that finds an organ-confined cancerous lesion.This radiological lesion is then confirmed as prostate cancer with GS = 3 + 3 in an ensuing biopsy.Because the prostate cancer risk is low, the patient enrolls in AS and periodic PSA tests are performed until the date of the next mpMRI scan in the AS protocol.This second imaging session does not reveal progression in the lesion, so the AS monitoring plan remains unchanged.However, the patient starts exhibiting a fast increase in PSA, which motivates an earlier imaging session before the originally prescribed date according to the AS protocol.This third mpMRI scan reveals radiological progression, which is further confirmed histopathologically as an upgrade to GS = 4 + 3 in an ensuing biopsy.At this point, the patient is offered a radical treatment for prostate cancer, which usually consists of surgery (i.e., radical prostatectomy) or radiotherapy (e.g., external beam radiotherapy, brachytherapy).B, Changes to the standard-of-care AS protocol after implementing the computational tumor forecasting pipeline presented in this study.The modified protocol is identical to the standard of care up to the second mpMRI.At this point, the longitudinal imaging data collected for the patient can be used to personalize our biomechanistic model of prostate cancer growth and obtain a computational forecast of prostate cancer growth over the patient's prostate anatomy.The forecast reveals progression toward high-risk prostate cancer and provides the time up to this event.This prediction enables optimizing the timing of the third mpMRI to confirm progression and proceed to treatment.Thus, our approach avoids PSA testing and biopsy after the second mpMRI, provides a personalized prediction of the patient's prostate cancer progression that enables an early detection of this event, and supports the decision and optimal timing to perform treatment.C, Main steps in our computational pipeline for prostate cancer forecasting during AS.In the cohort of this study (n = 7), all patients had three mpMRI scans.We first analyze the ability of our model to represent patient-specific prostate cancer growth after being informed by three mpMRI scans.Then, we also investigate the ability of our model to forecast prostate cancer growth when informed by only the first two mpMRI scans, and we use the third one to assess the predictive performance of the model.The first step of the computational pipeline is segmentation.We delineate the prostate and the tumor region of interest (ROI) on the longitudinal mpMRI data collected for the patient.After segmentation, the second and third mpMRI datasets and segments are co-registered with a nonrigid elastic method to the first one (registration transforms R 21 and R 31 , respectively).Next, we build a virtual model of the prostate anatomy, consisting of a 3D isogeometric (IGA) mesh and we project the registered tumor ROIs onto it.Then, we map the ADC values within each tumor ROI to tumor cell density values, which are subsequently used to guide model calibration and determine the personalized model parameters.Finally, we perform a patient-specific tumor forecast, including the prediction of tumor volume, tumor cell density map, and prostate cancer risk.

AACRJournals.org
Cancer Res Commun; 4(3) March 2024 619 interscan frequency of 1.4 (0.5, 2.6) years.For each patient, we define a relative time scale beginning at the time of the first mpMRI scan (t = 0).Although most biopsies were performed after an mpMRI scan, some patients also had biopsies before the first imaging session.Baseline GS at t = 0 was assigned according to the GS of the first biopsy performed within 3 months after the first mpMRI date.
If there was no biopsy within this period but the patient had a posterior biopsy (i.e., >3 months after the first mpMRI scan) and a preimaging biopsy, the first postimaging GS was used as the baseline measurement only if it was identical to the preimaging GS value.If the patient had no biopsies after the first mpMRI scan, the most recent GS measured before the first mpMRI was used as baseline.
All other postimaging GS measurements were used for temporally-informed clinical risk stratification during AS.This process enabled the definition of GS at t = 0 for all patients and resulted in a total of 16 eligible GS measurements with a median (range) of 3 (1,4) values per patient.These GS measurements were evenly distributed in two subgroups exhibiting a GS = 3 + 3 (n = 8) and a GS ≥ 3 + 4 (n = 8), which we denote as lower-risk and higher-risk prostate cancer.This division enabled a balanced subgroup definition and is consistent with the fact that prostate cancer with GS 3 + 3 is extensively recognized as eligible for AS, while there is debate regarding the inclusion in AS protocols of prostate cancer cases exhibiting Gleason grade 4 (2,12).

MRI Protocols
The mpMRI datasets were collected using 3T scanners for all subjects.A TIM Trio scanner (Siemens Healthineers USA) housed at the Research Imaging Institute of The University of Texas Health Science Center San Antonio (San Antonio, TX) was used to obtain a baseline scan with a cardiac coil used on the pelvis (37,38).Follow-up scans were performed following the standard-of-care using either a GE 750W (GE Healthcare USA) at the South Texas Veterans Healthcare System or a Siemens Skyra scanner (Siemens Healthineers USA) at The University of Texas Health Science Center San Antonio (San Antonio, TX).The mpMRI acquisition protocols included T 1 -weighted, T 2 -weighted, diffusion-weighted, and dynamic contrast-enhanced sequences (T1W, T2W, DW, and DCE, respectively).The main acquisition parameters were as follows: echo time <90 seconds, repetition time >3,000 ms, slice thickness <4.0 mm without gap, field of view of 160-220 mm, and b-values of approximately 50, 400, 800, and 1,400 s/mm 2 on all scanners.Apparent diffusion coefficient (ADC) maps were calculated from DW-MRI data using the standard monoexponential model (39).

Prostate Cancer Histopathology
All biopsies were performed using a transrectal ultrasound approach.After an mpMRI was performed, the regions of interest and prostate boundaries were marked in DynCAD software as part of the UroNav (Phillips Healthcare) biopsy system.During biopsy, once the prostate was scanned and registered with the MRI, three targeted biopsies were obtained from each region of interest and a standard 12-core prostate biopsy was also completed.Biopsy specimens underwent standard pathologic processing and were graded using the International Society of Urologic Pathology criteria and standard Gleason grade scoring.The same histopathologic assessment process was used for standard biopsies before the first mpMRI scan and for surgically-removed prostate specimens.

Preprocessing of Imaging Data
The integration of the mpMRI data from each patient within our biomechanistic model of prostate cancer growth required five preprocessing operations: (i) segmentation of the prostate in each T2W image, (ii) intrascan registration of the T2W map to the ADC map for each mpMRI scan, (iii) segmentation of the tumor on each ADC map, (iv) interscan registration of the longitudinal mpMRI data and the segmentations for each patient, and (v) conversion of intratumoral ADC values to tumor cell density estimates.This computational pipeline was designed following previous efforts developed to support mpMRIinformed personalized forecasting for brain (29) and breast cancers (21,28).For each mpMRI dataset, the segmentation of the prostate and the tumor, as well as the intrascan registration, were carried out using a combination of manual and semi-automatic tools available in 3DSlicer (40).Toward this end, the T2W and ADC maps were first resampled to isotropic resolution using linear interpolation.We next briefly describe the preprocessing operations in the order they were implemented.

Prostate Segmentation
The prostate gland was manually segmented in the resampled T2W images using anatomic landmarks as reference (3,4).To facilitate the ensuing intrascan registration, we also segmented the prostate in the ADC maps.

Intrascan Registration
To align the ADC map to the T2W map of each mpMRI dataset, we performed a rigid intrascan registration by leveraging the General Registration module of 3DSlicer and using the corresponding prostate segmentations for masking.Following this intrascan registration, we smoothed the T2W-based segmentations of the prostate using the Segmentation smoothing module in 3DSlicer.

Tumor Segmentation
The segmentation of the tumor consisted of a multistep, hierarchical procedure.
First, we manually drew a gross segmentation of the tumor over the registered ADC map based on the tumor location provided in the mpMRI reports, which was confirmed by an experienced urologist who specializes in prostate cancer and its detection using mpMRI (M.A. Liss).Second, we manually delineated a volume of healthy tissue in the same region of the prostate where the tumor is primarily located (i.e., peripheral zone or central gland) with approximately the same volume as the gross segmentation of the tumor.We then calculated the mean ADC over this volume as an estimate of ADC in healthy tissue (ADC h ) using the Segment Statistics module in 3DSlicer.Third, we used the threshold tool in 3DSlicer to delineate the tumor core within the gross segmentation of the tumor, which was defined as the isovolume of ADC values with less than 70% of ADC h .This value was selected after an analysis of five studies investigating the correlation of ADC and GS (41)(42)(43)(44)(45).These studies show that prostate cancer exhibits lower ADC than healthy prostatic tissue and that ADC progressively decreases as the GS increases.This trend seems to plateau for higher GS values, which aligns with a higher tumor cell density approaching tissue carrying capacity in those tumors (46)(47)(48).By fitting a hyperbolic tangent function to the average of the ADC ratios (ADC/ADC h ) obtained from the mean ADC h and mean GS-specific ADC values in the five studies (41)(42)(43)(44)(45), we determined that the 70% of ADC h was a sufficient threshold to identify mpMRI-visible tumors (see Fig. 2A).This fitting was performed using the Curve Fitting Toolbox in MATLAB (R2021b; The Mathworks) and defining a continuous GS scale ranging from 0 to 10, where values between 0 and 2 represent healthy tissue and pretumoral lesions (4).In addition, the choice of the hyperbolic tangent includes a slow decrease of ADC in pretumoral lesions and low GS tumors, which makes them practically indistinguishable from healthy tissue and aligns with the fact that mpMRI-observable prostate cancer generally exhibits GS ≥ 3 + 3 at diagnosis (3,14,16).Supplementary Data S1 provides further details about this hyperbolic tangent fit.Finally, the resulting tumor core was expanded to define the final tumor segmentation by leveraging the margin tool in 3DSlicer.This technique was applied within the gross tumor segmentation and using a  B).We use the linear mapping in Eq. ( 1), which introduces a negative proportionality between ADC and tumor cell density, as in other mpMRI-informed biomechanistic models of solid tumor growth (18,21,28,29).
margin of 2 to 4 mm (equivalent to a kernel from 3 × 3 × 3 to 5 × 5 × 5 voxels, depending on the tumor size), such that the tumor border lied over larger ADC values (e.g., 90% of ADC h ).Hence, this margin facilitated a smooth transition of ADC values from the tumor core to the neighboring healthy tissue that matches the natural solution of our biomechanistic model.

Interscan Registration
For each patient, the three mpMRI datasets collected during AS and their corresponding segmentations were co-registered to a common data frame, which was set as that of the first mpMRI scan (see Fig. 1C).This interscan registration was carried out using a biomechanically-constrained, deformable image registration algorithm to control for bulk deformations of the prostate anatomy between scan dates while preserving imaging features associated with longitudinal disease-related changes.The registration process consisted of a rigid registration of the T2W-based prostate segmentations followed by a deformable linearized iterative boundary reconstruction algorithm (49,50).This registration method filters out purely elastic deformation effects between imaging dates while preserving the effect of nonelastic imaging changes associated with tumor growth during AS.Further details are available in Supplementary Data S2.
The interscan registration was applied to the T2W data, ADC maps, and tumor segmentations of each patient.After the interscan registration, we updated the reference prostate segmentation (i.e., that defined in the T2W image of the first mpMRI scan) to prepare it for integration within our modeling framework by subtracting the urethral region.To do this, we used the eraser tool in 3DSlicer with a spherical geometry and 4-mm diameter.Anatomic landmarks in the T2W data were used for reference.

Conversion of Intratumoral ADC Values to Tumor Cell Density Estimates
Because ADC has been shown to decrease as tumor cell density increases in higher GS tumors (41)(42)(43)(44)(45)(46)(47)(48), the ADC values within each tumor segmentation were converted to normalized tumor cell density by leveraging a formulation that has been successfully employed to represent the inverse relationship between ADC and tumor cell density in other solid tumors (18,21,28,29), as follows: In Eq. ( 1), N(x, t ) is the normalized tumor cell density at position x and time t (i.e., 0 ≤ N(x, t ) ≤ 1), which is defined as the ratio of tumor cell density N(x, t ) to the tissue carrying capacity θ (i.e., the maximally admissible tumor cell density in the prostate, 0 ≤ N(x, t ) ≤ θ).In addition, ADC h is the ADC in the healthy tissue as defined above, ADC(x, t ) is the ADC at position x and time t within the tumor, and ADC min is the minimum ADC observable within a tumor.
For each ADC map, the latter was estimated from the lower horizontal asymptote of the hyperbolic tangent function in Fig. 2A as ADC min /ADC h = 0.25.
Figure 2B further provides the mapping through Eq. ( 1) of the hyperbolic tangent function describing the changes in ADC ratio with respect to GS, thereby yielding the changes in normalized tumor cell density with respect to GS within our modeling framework.Because ADC h and ADC min are fixed estimates, ADC values above ADC h and below ADC min were truncated to render a normalized tumor cell density of 0 and 1, respectively.This operation avoided unphysical negative tumor cell densities and values over the carrying capacity.Thus, according to Eq. ( 1) and the plots in Fig. 2, the segmentation of the tumor core using values below 70% ADC h corresponds to normalized tumor cell density values N ≥ 0.4, and, after the margin expansion, the border of the resulting tumor segmentation has values N 0.15.

Biomechanistic Model
We model prostate cancer growth during AS in terms of the spatiotemporal dynamics of the tumor cell density, N = N(x, t ), using a reaction-diffusion partial differential equation, as follows: Equation ( 2) is commonly known as the Fisher-Kolmogorov equation ( 51).Many patient-specific tumor forecasting studies have directly relied on this equation or leveraged it as a basis to build more sophisticated models (18,(21)(22)(23)(28)(29)(30)(31).The right-hand side of Eq. ( 2) describes the dynamics of the tumor cell density as a combination of two mechanisms: tumor cell mobility, which is represented by a diffusion process governed by the tumor cell diffusivity (D), and tumor cell net proliferation, which is modeled with a logistic term controlled by the net proliferation rate (ρ) and the tissue carrying capacity (θ).The tumor cell net proliferation rate encompasses the balance of tumor cell proliferation and death.In this work, we assume that the carrying capacity is constant for each patient, which enables posing our model in terms of the normalized tumor cell density ( N = N/θ) by dividing both sides of Eq. ( 2) by θ: Hence, once we calibrate the values of parameters D and ρ in Eq. ( 3 3) is that we want to leverage a parsimonious formulation of prostate cancer growth, such that it features a minimum set of personalized parameters that does not require extensive datasets, facilitates parameter identifiability, and avoids overfitting.Although θ may vary for each patient in this modeling framework, the normalization underlying Eq. ( 3) enables us to focus on relative changes of tumor cell density that match the changes in ADC observed in prostate cancer with increasing GS values in the literature (41)(42)(43)(44)(45). Thus, in our modeling framework θ can be patient-specific, but we do not need it to interrogate spatiotemporal prostate cancer dynamics.
Figure 3 illustrates how tumor cell mobility and net proliferation contribute to prostate cancer dynamics according to the models in Eqs. ( 2) and (3).As shown in Fig. 3, these models are posed within the patient's prostate geometry, which corresponds to the T2W-based segmentation from the first mpMRI scan.Furthermore, the model is defined over a temporal interval [0, T ], where T is a specific time horizon (e.g., for model calibration or tumor forecasting for each patient).To complete the model, we need to define boundary conditions over the prostate surface.Because the prostate cancer cases eligible for this study were organ-confined, we set zero-flux boundary conditions (i.e., ∇N • n = 0 or, equivalently, ∇ N • n = 0; where n is an outward unit vector orthogonal to the prostate boundary).In addition, the initial conditions for Eq.(3), N(x, 0), were always defined with the normalized tumor cell density map obtained from the first mpMRI scan.
To calibrate the model and obtain prostate cancer predictions for each patient, we performed computer simulations of our model by solving Eq. ( 3) numerically.These computer simulations constitute an in silico representation of personalized prostate cancer growth over the patient's prostate anatomy and the temporal interval [0, T ], which we can compare with mpMRI data at scan times.We leveraged isogeometric analysis (IGA), which is a high-fidelity gen-eralization of the finite element method that uses highly smooth functions to discretize the geometry and the solution of partial differential equations (18,52).In particular, we discretized Eq. ( 3) in space using a standard isogeometric Bubnov-Galerkin method that relied on a 3D C 1 quadratic non-uniform rational B-spline (NURBS) functional space (18,33,52).This method requires the construction of a patient-specific mesh of the prostate using the same function space.Toward this end, we used a parametric mapping method (17,32), which deformed a reference torus mesh to match the final prostate segmentation defined over the T2W image from the first mpMRI scan of each patient.The original torus and prostate meshes had a polar discretization with 32 elements in the circumferential directions and eight elements in the radial direction.Each prostate mesh was subsequently refined using standard knot insertion (52) for numerical accuracy in the computer simulations of our biomechanistic model.Hence, the final patient-specific prostate meshes had 128 elements in the circumferential directions and 32 elements in the radial direction.The tumor segmentations and corresponding normalized tumor cell density maps were then L  -projected over the prostate mesh for each patient (17,32,33).Furthermore, the time interval [0, T ] was discretized with a constant time step t = 1 day for all patients, which sufficed to capture the spatiotemporal dynamics observed in the patients of this study.We integrated in time using the generalized-α method (33,52).This iterative algorithm produces a series of nonlinear systems of equations in each time step, which we linearized with the Newton-Raphson method.The ensuing linear systems were solved utilizing the generalized minimal residual method.Further details on the implementation of these numerical methods to solve our biomechanistic model are available in refs.17, 18, 32, 33, 52.

Model Calibration
The personalized parameterization of our biomechanistic model consisted of finding the values of the parameter set (D, ρ) that minimize the relative squared difference between the model predictions and the corresponding imagingbased measurements of normalized tumor cell density at the times of the mpMRI scans used to inform the model.This calibration problem can be formulated as: where n s is the number of mpMRI scans used for model calibration, represents the 3D geometry of the patient's prostate, and Nd (x, t i ) and Nm (x, t i , D, ρ) denote the normalized tumor cell density measured from imaging data and calculated from the model with parameters D and ρ at the scan times t i (i = 1, . . ., n s ), respectively.Hence, Eq. ( 4) defines a nonlinear least-squares problem to find the patient-specific parameter values (D, ρ), which we solved by leveraging the Gauss-Newton method (53).The initial guess for the parameters was set as D = 5 • 10 Model performance during calibration and forecasting was assessed with a panel of global and local metrics.The global metrics were the tumor volume (V T ) and the total tumor cell volume (V N ), which are defined as: In Eqs. ( 5) and ( 6

Model-based Biomarkers of High-risk Prostate Cancer
We analyzed the personalized computer simulations of our model to find model-based biomarkers of high-risk prostate cancer.Toward this end, we first calculated the values of a panel of six candidate markers at the times of histopathologic assessment of the tumors in our patient cohort, and then compared their values between the subgroups of low-risk and high-risk prostate cancer (i.e., GS = 3 + 3 and GS ≥ 3 + 4, respectively).These candidate markers were calculated using the model simulations obtained with the personalized parameters from a global calibration study, in which the model was informed with the maximum amount of mpMRI data available for each patient in this study.The list of candidate markers was designed on the basis of previous tumor forecasting studies (22,23,(28)(29)(30)(31)(32)54) and our analysis of the personalized predictions in this work.The six candidate markers were: the prostate volume (V P ), the tumor volume (V T ), the total tumor cell volume (V N ), the mean normalized tumor cell density ( N), the total tumor index (N T ), and the mean proliferation activity of the tumor (A p ).The prostate volume and the last three markers were calculated as follows: The prostate volume and the total tumor index defined in Eqs. ( 7) and ( 9), respectively, were included to account for the organ size effect on tumor dynamics, because high-volume prostates enlarged by concomitant benign prostatic hyperplasia (BPH) have been suggested to exhibit a lower risk of high-grade prostate cancer (32,55).In addition, the mean proliferation activity of the tumor (A p ) was defined on the basis of the work presented in ref. 54.
Any candidate marker exhibiting a statistically significant difference between the low-risk and high-risk prostate cancer subgroups was considered an eligible model-based biomarker of high-risk prostate cancer.We analyzed the performance of these biomarkers by building univariate logistic regression classifiers and using ROC curve analysis.We further built bivariate logistic regression classifiers including a maximum of two markers from the candidate list, such that at least one was identified as a biomarker.We did not consider logistic regression classifiers with more than two input markers to avoid overfitting given the reduced size of the histopathologic measurements in the cohort (n = 16).All logistic regression classifiers were built using the function glmfit from the Statistics and Machine Learning Toolbox in MATLAB.The area under the ROC curve was calculated numerically with the trapezoidal rule and the optimal performance point was determined using the minimum distance to the upper left corner (i.e., where sensitivity and specificity are maximal).
Given the small size of the cohort in this work (n = 7), we did not partition it to define a patient subgroup to train the classifiers and another group to validate them.Nevertheless, to investigate their predictive performance, we analyzed them using the personalized computer simulations from the fitting-forecasting study.In particular, we analyzed whether the best performing classifier, which was trained with markers calculated with personalized simulations from the global calibration study, was able to anticipate the development of high-risk disease for each patient when the model was only informed by the first two mpMRI scans.

Statistical Methods
We leveraged Wilcoxon rank-sum and signed-rank tests to compare the local metrics of model performance calculated at the times of the second and third mpMRI scans in both the global calibration scenario and the fitting-forecasting study.These two statistical tests were also used to investigate differences in the model parameters as well as global and local metrics of model performance between the global calibration scenario and the fitting-forecasting study.In addition, the Wilcoxon rank-sum test was also used to identify differences between the candidate markers of high-risk prostate cancer between the low-risk and high-risk prostate cancer subgroups.The two types of statistical tests were performed using the functions ranksum and signrank from the Statistics and Machine Learning Toolbox in MATLAB.In the Results section, we specify when we use each type of test and whether it is two-tailed or one-tailed for each statistical analysis.The level of significance for all statistical tests was set to 5%.

Data Availability
The raw data leveraged in this study were generated at The University of Texas Health Science Center San Antonio (San Antonio, TX).Raw data are not publicly available under the IRB data usage agreement for this study, but the derived data that support the findings presented in this work are available from the corresponding author upon reasonable request.

The Biomechanistic Model Represents Prostate Cancer Growth During AS for Each Individual Patient
We first performed the global calibration study to assess the ability of our computational framework to represent untreated prostate cancer growth during AS.In this study, the biomechanistic model was informed by the three mpMRI scans available for each patient (n = 7).The reference prostate geometry from the first mpMRI scan had median (range) volume of 36.4 (18.5, 67.3) cc.The imaging-based segmentation of the tumors over the first, second, and third mpMRI datasets resulted in a tumor volume (V T ) of 0.12 (4.4 • 10

Personalized Computational Forecasts Predict Untreated Prostate Cancer Growth
To investigate the predictive performance of our computational framework, we carried out a fitting-forecasting study: for each patient, we first initialized our biomechanistic model with the data from the first mpMRI, we then  Similar results for the other 3 patients are provided in Supplementary Fig. S2.In addition, the unity plots in Fig. 7A and B compare the model predictions of the tumor volume (V T ) and the total tumor cell volume (V N ) to their corresponding imaging measurements.At the time of the second mpMRI scan (i.e., the calibration horizon), the PCC and CCC values for the model-data agreement of V T and V N were ≥0.99 in all cases.The PCC and the CCC for V T resulted in 0.97 and 0.96 at the third imaging timepoint (i.e., the forecasting horizon), respectively.The corresponding values of the PCC and the CCC for V N at the third mpMRI date were 0.93 in both cases.

Model-based Biomarkers Enable Early Identification of Progression to High-risk Prostate Cancer
To search for potential model-based biomarkers of high-risk prostate cancer, we investigated six quantities of interest that were calculated at the times of histopathologic assessment of each patient's tumor (n = 16) using the personalized model simulations from the global calibration study.The median (range) of the prostate volume (V P ), the tumor volume (V T ), the total tumor cell volume (V N ), the mean normalized tumor cell density ( N), the total tumor index (N T ),  8G shows the ROC curve for a univariate logistic classifier of high-risk prostate cancer that was constructed using the mean proliferation activity of the tumor (A p ).The AUC of this ROC curve is 0.77 and the optimal performance point operates at 75% sensitivity and specificity.In addition, Fig. 8G also shows the ROC curve for the best performing bivariate logistic classifier that could be built using the mean proliferation activity of the tumor (A p ) and one of the other markers, which resulted to be the total tumor index (N T ).The AUC for this ROC curve is 0.83 and the optimal performance point also operates at 75% sensitivity and specificity.
To investigate whether the bivariate logistic classifier would anticipate the detection of high-risk prostate cancer at the time of the second mpMRI scan (i.e., once the biomechanistic model is fully calibrated for each patient, and we can perform a personalized forecast), we tested it in the fitting-forecasting scenario.The results of this analysis for the 4 patients in Figs. 4 and 6 are shown in Fig. 9, while the corresponding results for the other 3 patients are reported in Supplementary Fig. S3.These figures further compare the performance of the classifier using the personalized model simulations from the global calibration scenario (i.e., those used for classifier training).The prostate cancer cases in Fig. 9A and D were consistently classified as low-risk and high-risk prostate cancer, respectively, in both the global calibration and the fitting-forecasting study, thereby matching the initial histopathologic assessment available for these tumors.Similar results were obtained for the patients with high-risk prostate cancer in Fig. 9B and Supplementary Fig. S3A.For the latter, our  (n = 8) and higher-risk (n = 8) prostate cancer, which were defined as tumors with GS = 3 + 3 and GS ≥ 3 + 4, respectively.These model-based biomarkers were calculated at the times were both a histopathologic assessment and imaging measurement are available for each patient in the cohort (n = 7).In particular, the model-based markers in A-F are: prostate volume (V P ), tumor volume, total tumor cell volume (V N ), mean normalized tumor cell density ( N), total tumor index (N T ), and mean proliferation activity of the tumor (A p ). Outliers are represented as hollow circles and an asterisk indicates significance under a one-sided Wilcoxon rank-sum test (P < 0.05).G, ROC curves for (i) the univariate logistic regression model constructed using the mean proliferation activity of the tumor (i.e., the only model-based marker that was significantly different between lower-risk prostate cancer and higher-risk prostate cancer), and (ii) the bivariate logistic regression model constructed using the mean proliferation activity of the tumor and the total tumor index (which is the combination of model-based markers that rendered the highest performance).The AUC of each ROC curve is reported within the plot, and the optimal performance point for both the univariate and bivariate logistic regression models operates at 75% sensitivity and specificity (bullet points on the ROC curves).and 6.In all panels, darker hues represent results from the global calibration scenario (GC, see Fig. 4), where the model is fit to the three mpMRI datasets from each patient, while lighter hues show results from the fitting-forecasting scenario (FF, see Fig. 6), where the model is only fit to the first two mpMRI datasets from each patient.Dotted gray vertical lines in the background indicate the times of the mpMRI scans for each patient.

AACRJournals.org
Solid lines correspond to quantities calculated from the model fit, while dashed lines correspond to values calculated from model forecasts.The model-based markers of interest are the total tumor index (N T , green curves, left vertical axis) and mean proliferation activity of the tumor (A p , pink curves, right vertical axis).The prostate cancer risk classifier was trained with the global calibration results (see Fig. 8), yielding an optimal performance threshold that separates lower risk prostate cancer (blue region) from higher-risk prostate cancer (red region).The prostate cancer risk at the times of histopathologic assessment of the tumors (i.e., biopsy, surgery) is represented as a bullet point, and the corresponding GS values are annotated over the prostate cancer risk trajectory from the global calibration scenario.In A, the patient exhibits a low-risk tumor during AS, which is correctly identified by our classifier in both computational scenarios.In B, our model consistently classifies the patient's tumor as higher-risk during the majority of AS in both computational scenarios.In C, the patient initially exhibits a lower-risk tumor that had progressed to higher-risk at the time of surgery (i.e., terminal GS value).In the fitting-forecasting scenario our model consistently predicts that the tumor is a higher-risk case, while, in the global calibration scenario, tumor progression is detected shortly after the second mpMRI scan.Importantly, in the fitting forecasting scenario, our approach identifies a higher-risk tumor 1,011 days earlier than standard practice (i.e., assessment of the surgical specimen).In D, our model consistently identifies the tumor as a higher-risk case in both computational scenarios.
modeling framework also anticipated the progression from GS 4 + 3 to GS 4 + 4 using the personalized predictions at the second imaging timepoint, which was 410 days (i.e., ∼1.1 years) earlier than the final histopathologic assessment at surgery.Importantly, although the early biopsies of the tumor in Fig. 9C report a low-risk case, the personalized biomechanistic model predictions (calculated at the second imaging timepoint) consistently classified this tumor as high risk and the global calibration (calculated at the third imaging timepoint) further confirmed prostate cancer progression right after the second scan.For this patient, our model anticipated the detection of prostate cancer progression at the time of the second mpMRI, which was 1,011 days (i.e., ∼2.8 years) before the final histopathologic assessment after surgery.Similar results are obtained for the prostate cancer case in Supplementary Fig. S3B, for which tumor progression was detected 677 days (i.e., ∼1.9 years) before surgery.Finally, the low-risk patient in Supplementary Fig. S3C was consistently identified as highrisk with the predictions of the fitting-forecasting study, but assimilation of the third mpMRI dataset in the global calibration scenario correctly classified the tumor as low-risk.

Discussion
The majority of patients with newly diagnosed prostate cancer are estimated to exhibit low or intermediate-risk prostate cancer at diagnosis (5,6), for which AS is a standard clinical option (2,3,12).To optimize patient monitoring during AS and rationalize the treatment of prostate cancer, it is of utmost importance to accurately identify the patients with indolent disease, who can continue to reap the benefits of AS and avoid treatment, as well as the patients with progressing prostate cancer, who may require immediate treatment.Nevertheless, current AS protocols largely rely on assessing prostate cancer status at standardized frequencies that are set upon the average dynamics of the disease observed in large clinical studies (2,12,13).This population-based, observational approach cannot anticipate prostate cancer progression and does not enable an optimal surveillance of the unique growth dynamics of each patient's tumor.To address these limitations of AS, we propose to use personalized computational forecasts of prostate cancer growth obtained with a spatiotemporal biomechanistic model informed by the longitudinal mpMRI data collected during AS for each individual patient.Hence, our forecasting technology ultimately aims at advancing AS from its current observational, population-based standard toward a predictive, patient-specific paradigm.
Our prostate cancer forecasts primarily consist of the prediction of spatiotemporal tumor growth in terms of normalized tumor cell density maps over the 3D anatomy of the patient's prostate.In the pilot study presented herein, we have shown that our approach can represent and predict prostate cancer growth according to global metrics of tumor volume and cellularity (see Figs. 5 and  7).This performance is analogous to that reported in prior tumor forecasting studies (18,28,29).Our results also show that our computational framework exhibits a promising potential to capture the evolving 3D morphology of the tumor within the patient's prostate (see Figs. 4 and 6).This information can be central to the precise planning of radical treatments, such as surgery or radiotherapy (17,18,22,23,29).In addition, computational predictions of the 3D tumor morphology within the patient's prostate can be leveraged to guide biopsies and, therefore, reduce the risk of underestimation of GS with respect to the histopathology analysis of the surgically-removed prostate (2,11).Furthermore, despite the simplicity of the biomechanistic model leveraged in this study, our results show a promising trend in the local agreement between imaging measurements and model calculations of the normalized tumor cell density maps (see Figs. 4-7), which is comparable to other tumor forecasting studies (18,23,28,29,31).
Our patient-specific forecasts further enable the prediction of prostate cancer risk by leveraging a set of biomarkers calculated from normalized tumor cell density maps obtained from the personalized model.In particular, in this pilot study we found that the mean proliferation activity of the tumor can be used as a promising biomarker of higher-risk prostate cancer, whose classifying performance improved in combination with the total tumor index (see Eqs. ( 9) and ( 10)).Tumor cell proliferation has been found to be a central mechanism driving tumor growth and treatment response in previous biomechanistic modeling studies, and model-based biomarkers of proliferation activity have been found to correlate with tumor progression and unfavorable pathologic outcomes (22,25,26,30,54,56).In the case of prostate cancer, previous clinical studies have also shown that a high proliferation activity has been associated with higher GS as well as with poorer prognosis, treatment outcomes, and survival (57,58).
Furthermore, the total tumor index combines three key variables in the analysis of prostate cancer growth in the clinical literature: cellularity (measured via the average normalized tumor cell density), tumor volume, and prostate volume (see Eq. ( 9)).While the first two have been directly associated with higher-risk prostate cancer (4,12,13,(41)(42)(43)(44)(45)(46)(47)(48) Another limitation is that our personalized forecasting approach does not account for uncertainty quantification and requires at least two mpMRI datasets to initialize and calibrate the model.For the pilot study presented in this work, we employed a deterministic approach following other previous studies pioneering the use of tumor forecasting in clinical settings (18,21,22,(28)(29)(30)(31).Hence, it suffices to work with three mpMRI datasets per patient.This amount of data enables the analysis of our tumor forecasting method in a global calibration scenario (in which we assess whether we can reproduce observed tumor growth and progression using all available mpMRI data) as well as in a fittingforecasting scenario (in which we assess whether we can predict tumor growth and progression at the time of the last mpMRI scan by informing the model with only the first two mpMRI datasets).We believe that the results obtained in these two scenarios provide a good initial assessment of our patient-specific tumor forecasting methods and motivate further development and validation over larger cohorts.Of note, despite the diverse interscan times during AS for each patient (which ranged from 0.5 to 2.6 years), we would like to remark that our methods show strong predictive capabilities regardless of the time between the mpMRI datasets used for calibration as well as of the time between the last mpMRI dataset used for the personalized calibration of our biomechanistic model horizon and the forecasting horizon (i.e., the times of the second and third mpMRI datasets for each patient in this study, respectively  (17), further studies are needed to investigate the joint calibration of the parameters of the equations governing the temporal dynamics of serum PSA and the spatiotemporal dynamics of tumor cell density, along with the analysis of their contributions to predict prostate cancer risk.In addition, genomic classifiers are gaining increased interest for risk stratification of newly diagnosed prostate cancer.These tests have been strongly and independently associated with adverse pathology, treatment failure, metastasis, cancer-specific and overall survival, and conversion from AS to definitive treatment (64)(65)(66).While their use is not standardized in AS protocols and they are seldom performed longitudinally, the results of these genomic classifiers could be combined with our model-based biomarkers to refine the prediction of prostate cancer progression during AS as well as to refine the construction of prior distributions of parameter values within a probabilistic forecasting framework.
Finally, the simplicity of the model is a central limitation of this study.In particular, the biomechanistic model used herein tends to homogenize the intratumoral spatial distribution of the normalized tumor cell density map.This effect is caused by the Fisher-Kolmogorov equation, which we used to describe untreated prostate cancer growth in our biomechanistic model (see Eqs.  and 6), as well as the relatively lower local model-data agreement with respect to previous studies employing more advanced models (18,23,28,29,31).In addition, the progressive model-driven spatial homogenization of the normalized tumor cell density maps precludes accurate predictions of local GS values, for example, using the hyperbolic tangent mapping built to support the imaging preprocessing pipeline (see Fig. 2).The analysis of these local GS maps would also require data on the spatial position of the prostate cancer-positive biopsy cores and surgical specimen regions, but this information was not available for this study.Nevertheless, the prediction of local GS maps would be especially important to precisely identify the regions of the tumor for biopsy (2,3,16), to adapt the treating fields during radiotherapy (22,23,29), and, in general, to provide richer predictions of prostate cancer progression during AS.Thus, future studies should investigate model extensions that incorporate mpMRI-informed spatial mechanisms that aim at retaining and evolving the intratumoral heterogeneity observed in imaging measurements of normalized tumor cell density maps.To achieve this goal, previous studies have leveraged heterogeneous and anisotropic parameterizations (18,21,23,28,29,31), mechanical constraints to tumor cell mobility and proliferation (18,21,28,29,31), and information on tumor-supporting vascularity (18,20,21,28,67).In particular, given that the prostate is a confined organ in the pelvic region and that both a tumor and potentially coexisting BPH may enlarge the prostate volume, the mechanical stress field within the organ can have a major impact on tumor dynamics (32).Thus, the extension of the model to include the mechanical effects on prostate cancer growth constitutes a promising avenue of research.Toward this end, it would also be necessary to obtain a segmentation of the central gland of the prostate in T2W images to define the mechanical stress due to BPH (32), which can be straightforwardly processed along the corresponding prostate segmentation within our computational pipeline.Another promising extension to refine the spatial heterogeneity of our prostate cancer forecasts would be the definition of a spatially-varying carrying capacity map (θ(x)) based on local measurements of vascularization obtained via dynamic contrast-enhanced MRI (67) and the local level of mechanical stress, which may limit tumor cell density via changes in prostatic tissue architecture caused by prostate cancer and BPH (4, 18, 32, 46-48, 55, 68, 69).
In the future, we believe that further developing our computational modeling framework addressing the aforementioned limitations will enable the construction of a digital twin (27,62) for each individual patient enrolling in AS.These technologies not only enable a seamless integration of patient data and personalized model forecasts, but also exploit them to define specific metrics that aid the treating physician in making decisions about the timing and type of further testing or treatments, thereby providing a robust and systematic framework for the optimal management of each patient's tumor.By further extending our biomechanistic model to describe the response to standard first-line treatments for prostate cancer after AS, the digital twin could also be trained to calculate therapeutic outcomes or expected survival, and, eventually, optimal regimens maximizing them while minimizing side effects (18,(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30)(31)(32)(33)62).Future technological developments in mpMRI can also contribute to new and more accurate imaging measurements of prostate cancer (3,70), which can contribute to better inform our model and enhance the reliability of our model predictions.Thus, we posit that the work presented in this study constitutes an important first step toward the construction of a digital twin enabling a personalized, predictive paradigm for the clinical management of prostate cancer.

Conclusion
We have presented a clinical-computational framework for the personalized forecast of untreated, newly diagnosed prostate cancer growth during AS constructed upon a spatiotemporal biomechanistic model informed by patientspecific longitudinal mpMRI data.This technology tracks and predicts the 3D development of the tumor within the patient's prostate geometry, while also enabling to anticipate the progression to higher risk-disease.Although further model development and investigation within larger cohorts are necessary, we posit that the personalized predictive technology presented in this pilot study can contribute to guide optimal clinical decision-making regarding the frequency and type of monitoring during AS, as well as the optimal time, modality, and planning of ensuing treatment.

FIGURE 2
FIGURE 2 Mapping ADC to N and GS.A, Nonlinear least-squares fit of a hyperbolic tangent model (blue solid line) of the mean values of ADC in the tumor across GS groups in previous studies (red dots), which were normalized with respect to the mean ADC in healthy tissue (ADC h )reported for each of them(41)(42)(43)(44)(45).We leverage a continuous, extended GS scale in which 0 < GS < 2 indicates healthy and pretumoral tissue and GS > 2 corresponds to prostate cancer (i.e., overlapping the standard discrete values of the GS used to assess histopathologic samples of prostate cancer; see refs.2, 4).Further details on the hyperbolic tangent model and the fitting method are reported in Supplementary Data S1.B, Mapping of the hyperbolic tangent model obtained for the normalized ADC values (ADC/ADC h , A) to normalized tumor cell density values (i.e., N, B).We use the linear mapping in Eq. (1),

FIGURE 3
FIGURE 3 Biomechanistic model of prostate cancer growth during AS.We describe the growth of newly diagnosed untreated prostate cancer in terms of the dynamics of the normalized tumor cell density, N(x, t ), which may vary between N = 0 far from the tumor to a maximum value of N = 1 inside the tumor.The dynamics of the normalized tumor cell density is governed by the two mechanisms shown in this figure: (i) tumor cell mobility, which expands the tumor and is represented with a diffusion process controlled by parameter D (i.e., the tumor cell mobility coefficient); and (ii) tumor cell net proliferation, which increases the local tumor cell density and is modeled as a logistic growth process controlled by parameter ρ (i.e., tumor cell net proliferation rate).Because newly diagnosed prostate cancer cases eligible for AS are organ-confined, we further assume that tumor cells cannot leave the patient's prostate by enforcing ∇ N • n = 0.
), T is the 3D tumor region where N ≥ Nth with Nth being a threshold value to systematically identify tumor tissue both in the imaging data and the model forecasts.Hence, this threshold facilitated the automatic delineation of a tumor region in the imaging data, for which the values over the tumor segmentation border were N 0.15, and from the model simulations, which provided a continuous N map within the range [0,1].For this purpose, we set Nth = 0.15 in all the calculations presented in this work.The global metrics V T and V N were calculated patientwise in each calibration scenario, and model-data agreement was assessed via the Pearson and concordance correlation coefficients over the cohort (PCC and CCC, respectively).The local metrics employed to analyze the performance of the model were: the Dice similarity coefficient (DSC) as well as the root mean squared error (RMSE), the local PCC, and the local CCC between the tumor cell density maps calculated from the imaging data and leveraging our biomechanistic model.These local metrics were calculated patientwise and we analyzed their distribution using boxplots.

FIGURE 4
FIGURE 4 Examples of personalized calculations of prostate cancer growth in the global calibration scenario.A-D, Prostate cancer growth at the dates of the second mpMRI scan during model calibration (t MRI,2 ; upper row in each panel) and the third mpMRI scan (t MRI,3 ; bottom row in each panel) in 4 patients.These results are illustrated using a 3D representation of the mpMRI-extracted prostate geometry of each patient including the imaging-measured and model-calculated tumor regions (red and blue volumes), along with an axial section of the prostate showing the normalized tumor cell density map obtained from the mpMRI data and using the personalized model (i.e., Eq. (3)).

Figure 4
Figure 4 illustrates the 3D imaging measurements and calculations of the tumor geometry and the normalized tumor cell density map obtained with the personalized model for 4 patients at the second and third imaging timepoints.The corresponding results for the other 3 patients are provided in Supplementary Fig. S1.In addition, Fig.5A and Bprovide unity plots comparing the imaging measurements and the model calculations of the tumor volume (V T ) and the total tumor cell volume (V N ).Considering the second mpMRI scan for each patient, the PCC and CCC for V T were 0.99 and 0.97, respectively.The corresponding PCC and CCC for V N were 0.99 and 0.98.At the time of the third mpMRI scan, the PCCs and CCCs of both V T and V N were ≥0.99 in all

FIGURE 5
FIGURE 5 Summary of metrics assessing the performance of the biomechanistic model during global calibration.A and B, Unity plots to assess the agreement between the model estimation and mpMRI measurement of the tumor volume (V T ) and the total tumor cell volume (V N ) across the patient cohort (n = 7).These global tumor metrics were calculated at the dates of the second and third mpMRI scans (blue and green points).The unity plots in A and B also report the corresponding PCC and CCC for the global tumor metrics at the dates of the second and third mpMRI scans.C-F, Distributions of four local metrics assessing the agreement between the tumor cell density map calculated with our biomechanistic model and extracted from mpMRI data across the patient cohort in the global calibration scenario: Dice similarity coefficient (DSC), root mean squared error (RMSE), Pearson correlation coefficient (PCC), and concordance correlation coefficient (CCC).These local metrics were calculated at the dates of the second and third mpMRI scans (blue and green boxplots).

Figure 6
Figure6presents the model-predicted tumor geometry and normalized tumor cell density map obtained in the fitting-forecasting study at the times of the second and third mpMRI scans for the same 4 patients shown in Fig.4in the global calibration study, along with the corresponding imaging-based measurements.
Figure7C-E also present the distributions of the local metrics quantifying the pointwise agreement between the imaging measurements and model forecasts of the normalized tumor cell density maps ( N(x, t )) over the patient cohort (n = 7).At the second imaging timepoint, we obtained a DSC of 0.82 (0.67, 0.85), an RMSE of 0.022 (0.008, 0.068), a local PCC of 0.55 (0.19, 0.67), and a local CCC of 0.55 (0.18, 0.66).The corresponding values of the DSC, RMSE, local PCC, and local CCC at the forecasting horizon were 0.76 (0.65, 0.85), 0.048 (0.016, 0.137), 0.37 (0.26, 0.61), and 0.33 (0.24, 0.60), respectively.Two-sided Wilcoxon rank-sum testing resulted in no significant differences (P > 0.05) between the values of the local metrics at the times of the second and third mpMRI scans, although the DSC values were identified as significantly higher at the calibration horizon with respect to the forecasting timepoint under a one-sided Wilcoxon rank-sum test (P = 0.049).According to two-sided Wilcoxon signed-rank tests, model forecasts at the third imaging timepoint exhibited a significantly higher RMSE (P = 0.016) and significantly lower DSC, local PCC, and local CCC (P = 0.016, 0.047, and 0.047, respectively).No significant differences were detected between the values of the model parameters (D, ρ) obtained in the global calibration and the fitting-forecasting study under both Wilcoxon rank-sum and signed-rank testing (P > 0.05).Nevertheless, we observed a tendency toward lower tumor cell diffusivity coefficients (D) and higher net tumor cell proliferation rates (ρ) in the fittingforecasting study for each patient under one-sided Wilcoxon signed-rank tests (P = 0.055).In addition, no Wilcoxon rank-sum test identified significant differences between the local model-data agreement metrics from the global calibration and the fitting-forecasting studies at either the second or third imaging timepoints (P > 0.05).However, two-sided Wilcoxon signed-rank testing resulted in significantly lower RMSE values at the third mpMRI date in the global calibration study (P = 0.016) and at the second imaging timepoint in the fitting-forecasting study (P = 0.016).Furthermore, the DSC at the second scan timepoint was significantly higher in the fitting-forecasting study under a one-sided Wilcoxon signed-rank test (P = 0.039).

FIGURE 6
FIGURE 6 Examples of personalized calculations of prostate cancer growth in the fitting-forecasting scenario.A-D, Prostate cancer growth at the date of the second mpMRI scan during model calibration (t MRI,2 ; upper row in each panel) and the ensuing prediction of tumor growth at the date of the third mpMRI scan (t MRI,3 ; bottom row in each panel) in 4 patients.The patients shown in this figure are the same considered in Fig. 4.These results are illustrated using a 3D representation of the mpMRI-extracted prostate geometry of each patient including the imaging-measured and modelcalculated tumor regions (red and blue volumes), along with an axial section of the prostate showing the normalized tumor cell density map obtained from the mpMRI data and using the personalized model (i.e., Eq. (3)).

FIGURE 7
FIGURE 7 Summary of metrics assessing the performance of the biomechanistic model in the fitting-forecasting scenario.A and B, Unity plots to assess the agreement between the model estimation and mpMRI measurement of the tumor volume (V T ) and the total tumor cell volume (V N ) across the patient cohort (n = 7).These global tumor metrics were calculated at the dates of the second and third mpMRI scans (blue and green points).The unity plots in A and B further report the corresponding PCC and CCC for the global tumor metrics at the dates of the second and third mpMRI scans.C-F, Distributions of four local metrics assessing the agreement between the tumor cell density map calculated with our biomechanistic model and extracted from mpMRI data across the patient cohort in the global calibration scenario: Dice similarity coefficient (DSC), root mean squared error (RMSE), Pearson correlation coefficient (PCC), and concordance correlation coefficient (CCC).These local metrics were calculated at the dates of the second and third mpMRI scans (blue and green boxplots).

FIGURE 8
FIGURE 8 Potential model-based biomarkers of higher-risk prostate cancer.A-F, Distribution of six potential model-based biomarkers in lower-risk

FIGURE 9
FIGURE 9 Patient-specific forecasts of prostate cancer risk.A-D, Time trajectories of the model-based markers involved in the calculation of our prostate cancer risk classifier (left) as well as the trajectory of the latter (right) for 4 patients.The patients shown in this figure are the same considered in Figs.4 and 6.In all panels, darker hues represent results from the global calibration scenario (GC, see Fig.4), where the model is fit to the three mpMRI datasets from each patient, while lighter hues show results from the fitting-forecasting scenario (FF, see Fig.6), where the model is only fit to the first two mpMRI datasets from each patient.Dotted gray vertical lines in the background indicate the times of the mpMRI scans for each patient.
(2) and(3)).The natural solution of this equation is a travelling wavefront (51) that volumetrically expands a spatial region with nonzero tumor cell density to occupy other neighboring healthy regions, such that normalized tumor cell density progressively increases in the central part of the tumor and decreases along the travelling wavefront.Hence, although the initial conditions defined upon the normalized tumor cell density map extracted from the first mpMRI of each patient introduce information about the internal heterogeneity of the tumor architecture, the ensuing model simulation will progressively adapt the internal spatial distribution of normalized tumor cell density to match the natural solution of the Fisher-Kolmogorov equation.This effect contributes to explain the comparatively poorer performance of the personalized model at the third imaging timepoint with respect to the second one in both the global calibration and fitting-forecasting studies (see Figs.4 the mean proliferation activity of the tumor (A p ) in the low-risk prostate cancer group (i.e., GS = 3 + 3; n = 8) were 37.7 (21.8, 39.5) cc, 0.17 (0.01, 5.09) cc, 0.07 (1.5 • 10 −3 , 2.12) cc, 0.39 (0.26, 0.52), 2.31 • 10 −3 (3.82 • 10 −5 , 5.36 • 10 −2 ), and high-risk prostate cancer subgroups.Only the mean proliferation activity of the tumor (A p ) was significantly larger in high-risk prostate cancer according to a one-sided Wilcoxon rank-sum test (P = 0.041).Figure (27,(59)(60)(61) been suggested to have an inverse correlation (i.e., such that larger prostates tend to harbor tumors with more pathologically favorable features; see refs.32,55).Thus, the model-based quantities selected in the current study to construct a predictive classifier of high-risk prostate cancer agree with previous analyses of the driving mechanisms of tumor growth in modeling studies, as well as observations of prostate cancer biology and clinical progression in the literature.These results along with the promising performance of the bivariate classifier to predict progression at the time of the second mpMRI scan suggest that our forecasts of prostate cancer risk could contribute to guide fundamental clinical decisions in AS, such as the frequency of monitoring tests and the best timing to direct a patient with progressing prostate cancer to definitive treatment.For example, should the prostate cancer risk prediction reveal that a tumor will be at low risk at the 1-to 3-year horizon, then the interscan time can be extended to monitor the tumor.Conversely, if the forecast of prostate cancer risk reveals a trend toward progression, the subsequent scans could be performed earlier to confirm this pathologic event and promptly treat the tumor.Hence, only the longitudinal data from each individual patient is needed to constrain the model and make personalized forecasts, whose validation only depends on the subsequent mpMRI data collected from the patient during the course of AS.While a large cohort may provide moreAACRJournals.orgCancerResCommun;4(3)March2024 629 information of our biomechanistic model's predictive performance over more diverse prostate cancer cases, the calculation and validation of additional forecasts of spatiotemporal growth of prostate cancer would still be patient-specific (i.e., our biomechanistic model does not necessitate a cohort of patients for training, contrary to pure data-driven approaches(27,(59)(60)(61)).Conversely, our predictions of prostate cancer progression rely on a data-driven approach (i.e., logistic classifiers) and, given the small size of the patient cohort used in our study, we did not define separate training and validation cohorts.Despite this limitation, the promising preliminary results on prostate cancer progression presented in this work now justify performing subsequent studies in large patient populations featuring different training and validation cohorts to assess the performance, reproducibility, and reliability of our personalized forecasts and model-based biomarkers.Second, we used both biopsy and surgicallybased measurements of GS to train our high-risk prostate cancer classifiers.However, prostate cancer biopsies may underestimate GS mainly due to sampling bias and the incomplete understanding of imaging-histopathology data correlation(3, 11-13, 16, 48).Thus, a more robust training of our classifiers could be achieved by exclusively leveraging GS measurements on surgically removed prostate specimens paired to recent presurgery mpMRI scans.
(17,32)3))less, to maximize the utility of our predictive technology during AS, it would be necessary to enable tumor forecasting from the first imaging timepoint and then update the model as more data become available.Indeed, the comparison of the results of the global calibration study (three mpMRI scans per patient) and those from the fitting-forecasting study (two imaging datasets per patient) In such a probabilistic framework, the model provides a collection of tumor forecasts for multiple samples of the distribution of the model parameters based on existing data for each patient and prior knowledge over the population, thereby enabling predictions from baseline.Then, these forecasts can be analyzed using standard statistical metrics (e.g., median and range, mean and 95% confidence interval) and risk-based metrics (e.g., quantiles and superquantiles) that enable uncertainty quantification.In particular, the latter can inform about the risk of tumor cell density or other tumor biomarker reaching thresholds suggesting tumor progression(62).Furthermore, as more longitudinal data are collected for the patient, the parameter distributions can be updated, and this additional information can be transferred to update the distributions of tumor forecasts and biomarkers via Bayesian data assimilation(27,62,63).Because reproducibility and reliability are substantial barriers for the clinical adoption of new imaging biomarkers, uncertainty quantification could also contribute to the clinical translation of the model-based biomarkers proposed herein under variations in image acquisition as well as segmentation and registration errors.Our model is only informed by mpMRI data, which constitutes another limitation of this work.Despite the wealth of spatiotemporal information provided by this imaging modality, patients with prostate cancer in AS are also monitored using PSA tests, which are often performed more frequently than mpMRI scans(2,12,13).Thus, combining the spatiotemporal model of tumor cell dynamics used here (i.e., Eqs.(2) or (3)) with an equation describing PSA dynamics could facilitate a more frequent update of the patient-specific parameters and include PSA dynamics in the construction of the prostate cancer risk classifier(17,32).
tification of input data and model outcomes (e.g., parameters, model-based biomarkers), the update of the distributions of model parameters and predictions with new data collected during the course of AS, and tumor forecasting from the first imaging session.The reason for not including serum PSA in the current model is 2-fold.First, we wanted to build a parsimonious model that could reliably reproduce and predict spatiotemporal prostate cancer growth, for which mpMRI data suffice, and progression, for which GS data are necessary.Second, serum PSA is a scalar value and tumor cell density is a spatial quantity.Although we have already proposed a method to relate these two concepts